Click on tabs to display additional information.
# Load required libraries for statistical analysis and table creation
library(dplyr)
library(ggplot2)
library(gt)
library(broom)
library(car)
library(emmeans)
library(multcomp)
library(MASS) # For negative binomial regression
library(rstatix)
library(coin)
library(phyloseq)
library(microViz)
library(tidyverse)
# Define treatment order and color palette
treatment_order <- c(
"A- T- P-", # Control
"A- T- P+", # Parasite
"A+ T- P-", # Antibiotics
"A+ T- P+", # Antibiotics_Parasite
"A- T+ P-", # Temperature
"A- T+ P+", # Temperature_Parasite
"A+ T+ P-", # Antibiotics_Temperature
"A+ T+ P+" # Antibiotics_Temperature_Parasite
)
# Custom color palette matching treatment order
treatment_colors <- c(
"#1B9E77", # A- T- P- (Control)
"#D95F02", # A- T- P+ (Parasite)
"#7570B3", # A+ T- P- (Antibiotics)
"#E7298A", # A+ T- P+ (Antibiotics_Parasite)
"#66A61E", # A- T+ P- (Temperature)
"#E6AB02", # A- T+ P+ (Temperature_Parasite)
"#A6761D", # A+ T+ P- (Antibiotics_Temperature)
"#666666" # A+ T+ P+ (Antibiotics_Temperature_Parasite)
)
# Create named vector for color scale
treatment_color_scale <- setNames(treatment_colors, treatment_order)
# Function to extract sample data as dataframe from phyloseq object
samdatAsDataframe <- function(ps) {
samdat <- phyloseq::sample_data(ps)
df <- data.frame(samdat, check.names = FALSE, stringsAsFactors = FALSE)
return(df)
}
# Function to rename variables in phyloseq object
ps_rename <- function(ps, ...) {
ps <- microViz::ps_get(ps)
df <- samdatAsDataframe(ps)
df <- dplyr::rename(.data = df, ...)
phyloseq::sample_data(ps) <- df
return(ps)
}
ps.unclean <- readRDS('/Users/michaelsieler/Dropbox/Mac (2)/Documents/Sharpton_Lab/Projects_Repository/Rules_of_Life/major-experiment-2023/Data/Robjects/pseq_uncleaned_05052025.rds')
# Create a dataframe for mortality using the "unclean" phyloseq object
# The unclean phyloseq object is prior to removing samples during post-DADA2 processing (e.g., insufficient reads/sample)
data.mortality <-
ps.unclean %>%
## Update Metadata
ps_rename(Time = Timepoint) %>%
microViz::ps_mutate(
Treatment = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "A- T- P-",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "A- T- P+",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "A+ T- P-",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "A+ T- P+",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "A- T+ P-",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "A- T+ P+",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "A+ T+ P-",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "A+ T+ P+",
TRUE ~ "Unknown"
), .after = "Pathogen"
) %>%
microViz::ps_mutate(Sample = fecal.sample.number, .before = 1) %>%
microViz::ps_mutate(Sample = gsub("^f", "", Sample)) %>%
microViz::ps_filter(Treatment != "Unknown") %>%
microViz::ps_mutate(
History = case_when(
Antibiotics + Temperature == 0 ~ 0,
Antibiotics + Temperature == 1 ~ 1,
Antibiotics + Temperature == 2 ~ 2,
), .after = "Treatment"
) %>%
## Additional metadata updates, factorizing metadata
microViz::ps_mutate(
# Create treatment code
treatment_code = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 0 ~ "Aneg_Tneg_Pneg",
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Aneg_Tneg_Ppos",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Apos_Tneg_Pneg",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Apos_Tneg_Ppos",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Aneg_Tpos_Pneg",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Aneg_Tpos_Ppos",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Apos_Tpos_Pneg",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Apos_Tpos_Ppos"
),
# Create treatment group factor
treatment_group = case_when(
Antibiotics == 0 & Temperature == 0 & Pathogen == 1 ~ "Parasite",
Antibiotics == 1 & Temperature == 0 & Pathogen == 0 ~ "Antibiotics",
Antibiotics == 1 & Temperature == 0 & Pathogen == 1 ~ "Antibiotics_Parasite",
Antibiotics == 0 & Temperature == 1 & Pathogen == 0 ~ "Temperature",
Antibiotics == 0 & Temperature == 1 & Pathogen == 1 ~ "Temperature_Parasite",
Antibiotics == 1 & Temperature == 1 & Pathogen == 0 ~ "Antibiotics_Temperature",
Antibiotics == 1 & Temperature == 1 & Pathogen == 1 ~ "Antibiotics_Temperature_Parasite",
TRUE ~ "Control"
),
# Convert to factor with appropriate levels
treatment_group = factor(treatment_group,
levels = c("Control", "Parasite",
"Antibiotics", "Antibiotics_Parasite",
"Temperature", "Temperature_Parasite",
"Antibiotics_Temperature", "Antibiotics_Temperature_Parasite")
),
treatment_code = factor(treatment_code, levels = treatment_order),
# Create time point factor
time_point = factor(Time, levels = c(0, 14, 18, 25, 29, 60)),
# Create pathogen status factor
pathogen_status = factor(ifelse(Pathogen == 1, "Exposed", "Unexposed"),
levels = c("Unexposed", "Exposed")),
# Create sex factor
sex = factor(Sex, levels = c("M", "F"))
) %>%
microViz::ps_mutate(Treatment = factor(Treatment, levels = treatment_order)) %>%
microViz::ps_mutate(Exp_Type = case_when(
Treatment %in% c("A- T- P-", "A- T- P+") ~ "No prior stressor(s)",
Treatment %in% c("A+ T- P-", "A+ T- P+") ~ "Antibiotics",
Treatment %in% c("A- T+ P-", "A- T+ P+") ~ "Temperature",
Treatment %in% c("A+ T+ P-", "A+ T+ P+") ~ "Combined",
)) %>%
microViz::ps_mutate(Exp_Type = factor(Exp_Type, levels = c("No prior stressor(s)", "Antibiotics", "Temperature", "Combined"))) %>%
microViz::samdat_tbl()
data.infection <- data.mortality %>%
dplyr::filter(Pathogen == 1) %>%
dplyr::group_by(Treatment) %>%
dplyr::summarise(
n_fish = n(),
n_infected = sum(Total.Worm.Count > 0),
mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
percent_infected = round((n_infected / n_fish) * 100, 1),
.groups = "drop"
) %>%
dplyr::mutate(
# Convert Treatment to factor with specific order
Treatment = factor(Treatment, levels = treatment_order)
)
data.infection_tank <- data.mortality %>%
dplyr::filter(Pathogen == 1) %>%
dplyr::group_by(Treatment, Tank.ID) %>%
dplyr::summarise(
n_fish = n(),
n_infected = sum(Total.Worm.Count > 0),
mean_worm_burden = mean(Total.Worm.Count, na.rm = TRUE),
percent_infected = round((n_infected / n_fish) * 100, 1),
.groups = "drop"
) %>%
dplyr::mutate(
# Convert Treatment to factor with specific order
Treatment = factor(Treatment, levels = treatment_order)
)
# data.mortality
# data.infection
# data.infection_tank
Click on tabs to display tables. Scroll to see additional rows.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# Prepare data for statistical analysis - percent mortality rates
mortality_stats_data <- data.mortality %>%
dplyr::group_by(Treatment) %>%
dplyr::count() %>%
dplyr::mutate(
total_per_group = 90,
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1)
)
# Chi-square test for independence between treatment and mortality
set.seed(42)
chi_square_test <- chisq.test(
matrix(c(mortality_stats_data$dead, mortality_stats_data$n),
ncol = 2, byrow = FALSE)
)
# Create mortality rate summary table for display
mortality_rate_table <- mortality_stats_data %>%
dplyr::select(Treatment, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Summary - All Treatments",
subtitle = "Percent mortality and survival rates by treatment"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Chi-square test results
test_stat_table <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(chi_square_test$statistic, 4),
chi_square_test$parameter,
round(chi_square_test$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - All Treatments Mortality",
subtitle = "Test for independence between treatment and mortality"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
pairwise_results <- list()
treatments <- mortality_stats_data$Treatment
for(i in 1:(length(treatments)-1)) {
for(j in (i+1):length(treatments)) {
trt1 <- treatments[i]
trt2 <- treatments[j]
# Get data for these two treatments
data1 <- mortality_stats_data %>% dplyr::filter(Treatment == trt1)
data2 <- mortality_stats_data %>% dplyr::filter(Treatment == trt2)
# Create 2x2 contingency table
cont_table <- matrix(
c(data1$dead, data1$n, data2$dead, data2$n),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(trt1, trt2)
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
Comparison = paste(trt1, "vs", trt2),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Combine pairwise results
pairwise_summary <- do.call(rbind, pairwise_results) %>%
dplyr::arrange(desc(Significant), P_value) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - All Treatments Mortality",
subtitle = "Post-hoc comparisons between treatment pairs"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Show only significant pairwise comparisons
pairwise_summary__SigOnly <- do.call(rbind, pairwise_results) %>%
dplyr::arrange(desc(Significant), P_value) %>%
dplyr::filter(Significant == "Yes") %>%
gt::gt() %>%
gt::tab_header(
title = "Significant Pairwise Chi-Square Tests - All Treatments Mortality",
subtitle = "Post-hoc comparisons between treatment pairs (significant only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Rate Summary - All Treatments | ||||
| Percent mortality and survival rates by treatment | ||||
| dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|
| A- T- P- | ||||
| 9 | 81 | 10.0 | 90 | 90.0 |
| A- T- P+ | ||||
| 8 | 82 | 8.9 | 90 | 91.1 |
| A+ T- P- | ||||
| 14 | 76 | 15.6 | 90 | 84.4 |
| A+ T- P+ | ||||
| 10 | 80 | 11.1 | 90 | 88.9 |
| A- T+ P- | ||||
| 19 | 71 | 21.1 | 90 | 78.9 |
| A- T+ P+ | ||||
| 13 | 77 | 14.4 | 90 | 85.6 |
| A+ T+ P- | ||||
| 30 | 60 | 33.3 | 90 | 66.7 |
| A+ T+ P+ | ||||
| 20 | 70 | 22.2 | 90 | 77.8 |
| Chi-Square Test Results - All Treatments Mortality | |
| Test for independence between treatment and mortality | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 29.797800 |
| Degrees of freedom | 7.000000 |
| P-value | 0.000103 |
| Pairwise Chi-Square Tests - All Treatments Mortality | |||
| Post-hoc comparisons between treatment pairs | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| A- T- P+ vs A+ T+ P- | 14.7109 | 0.000125 | Yes |
| A- T- P- vs A+ T+ P- | 13.0933 | 0.000296 | Yes |
| A+ T- P+ vs A+ T+ P- | 11.6036 | 0.000658 | Yes |
| A- T+ P+ vs A+ T+ P- | 7.8221 | 0.005161 | Yes |
| A+ T- P- vs A+ T+ P- | 6.7680 | 0.009280 | Yes |
| A- T- P+ vs A+ T+ P+ | 5.1175 | 0.023686 | Yes |
| A- T- P+ vs A- T+ P- | 4.3573 | 0.036851 | Yes |
| A- T- P- vs A+ T+ P+ | 4.1105 | 0.042617 | Yes |
| A- T- P- vs A- T+ P- | 3.4258 | 0.064187 | No |
| A+ T- P+ vs A+ T+ P+ | 3.2400 | 0.071861 | No |
| A- T+ P- vs A+ T+ P- | 2.8042 | 0.094019 | No |
| A+ T- P+ vs A- T+ P- | 2.6307 | 0.104813 | No |
| A+ T+ P- vs A+ T+ P+ | 2.2431 | 0.134214 | No |
| A- T+ P+ vs A+ T+ P+ | 1.3358 | 0.247775 | No |
| A- T- P+ vs A+ T- P- | 1.2946 | 0.255204 | No |
| A- T+ P- vs A- T+ P+ | 0.9502 | 0.329676 | No |
| A+ T- P- vs A+ T+ P+ | 0.9065 | 0.341038 | No |
| A- T- P+ vs A- T+ P+ | 0.8625 | 0.353031 | No |
| A- T- P- vs A+ T- P- | 0.7976 | 0.371823 | No |
| A+ T- P- vs A- T+ P- | 0.5937 | 0.440995 | No |
| A- T- P- vs A- T+ P+ | 0.4661 | 0.494809 | No |
| A+ T- P- vs A+ T- P+ | 0.4327 | 0.510671 | No |
| A+ T- P+ vs A- T+ P+ | 0.1994 | 0.655213 | No |
| A- T- P+ vs A+ T- P+ | 0.0617 | 0.803785 | No |
| A- T- P- vs A- T- P+ | 0.0000 | 1.000000 | No |
| A- T- P- vs A+ T- P+ | 0.0000 | 1.000000 | No |
| A+ T- P- vs A- T+ P+ | 0.0000 | 1.000000 | No |
| A- T+ P- vs A+ T+ P+ | 0.0000 | 1.000000 | No |
| Significant Pairwise Chi-Square Tests - All Treatments Mortality | |||
| Post-hoc comparisons between treatment pairs (significant only) | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| A- T- P+ vs A+ T+ P- | 14.7109 | 0.000125 | Yes |
| A- T- P- vs A+ T+ P- | 13.0933 | 0.000296 | Yes |
| A+ T- P+ vs A+ T+ P- | 11.6036 | 0.000658 | Yes |
| A- T+ P+ vs A+ T+ P- | 7.8221 | 0.005161 | Yes |
| A+ T- P- vs A+ T+ P- | 6.7680 | 0.009280 | Yes |
| A- T- P+ vs A+ T+ P+ | 5.1175 | 0.023686 | Yes |
| A- T- P+ vs A- T+ P- | 4.3573 | 0.036851 | Yes |
| A- T- P- vs A+ T+ P+ | 4.1105 | 0.042617 | Yes |
Click on tabs to display tables. Scroll to see additional rows.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# Prepare data for statistical analysis - parasite exposure comparison (percent mortality rates)
parasite_mortality_data <- data.mortality %>%
dplyr::filter(Treatment %in% c("A- T- P-", "A- T- P+")) %>%
dplyr::group_by(Treatment) %>%
dplyr::count() %>%
dplyr::mutate(
total_per_group = 90,
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1)
)
# Fisher's exact test for 2x2 contingency table (more appropriate for small samples)
set.seed(42)
fisher_test <- fisher.test(
matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n),
ncol = 2, byrow = FALSE)
)
# Chi-square test as alternative
chi_square_parasite <- chisq.test(
matrix(c(parasite_mortality_data$dead, parasite_mortality_data$n),
ncol = 2, byrow = FALSE)
)
# Create mortality rate summary table for display
parasite_mortality_rate_table <- parasite_mortality_data %>%
dplyr::select(Treatment, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Summary - Parasite Exposure",
subtitle = "Percent mortality and survival rates: A-T-P- vs A-T-P+"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Fisher's exact test results
fisher_results <- data.frame(
Statistic = c("Fisher's exact test p-value", "Odds ratio", "95% CI lower", "95% CI upper"),
Value = c(
round(fisher_test$p.value, 6),
round(fisher_test$estimate, 4),
round(fisher_test$conf.int[1], 4),
round(fisher_test$conf.int[2], 4)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Fisher's Exact Test Results - Parasite Exposure Mortality",
subtitle = "Comparison of A-T-P- vs A-T-P+ mortality rates"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Chi-square test results for comparison
chi_square_parasite_results <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(chi_square_parasite$statistic, 4),
chi_square_parasite$parameter,
round(chi_square_parasite$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - Parasite Exposure Mortality",
subtitle = "Alternative test for A-T-P- vs A-T-P+ mortality rates"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Additional analysis: Compare mortality rates between parasite-exposed and unexposed
# Calculate difference in mortality rates
parasite_mortality_diff <- parasite_mortality_data %>%
dplyr::arrange(Treatment) %>%
dplyr::mutate(
Pathogen_Status = dplyr::case_when(
stringr::str_detect(Treatment, "P\\+") ~ "Exposed",
stringr::str_detect(Treatment, "P\\-") ~ "Unexposed",
TRUE ~ "Unknown"
)
) %>%
dplyr::select(Pathogen_Status, percent_mortality)
## Adding missing grouping variables: `Treatment`
# Calculate the difference in mortality rates
mortality_rate_diff <- parasite_mortality_diff$percent_mortality[2] - parasite_mortality_diff$percent_mortality[1]
# Create summary of mortality rate difference
parasite_mortality_diff_summary <- data.frame(
Comparison = "A-T-P+ vs A-T-P-",
Mortality_Rate_Difference = round(mortality_rate_diff, 2),
Unexposed_Rate = round(parasite_mortality_diff$percent_mortality[1], 2),
Exposed_Rate = round(parasite_mortality_diff$percent_mortality[2], 2),
Percent_Increase = round((mortality_rate_diff / parasite_mortality_diff$percent_mortality[1]) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Difference - Parasite Exposure",
subtitle = "Comparison of mortality rates between parasite-exposed and unexposed groups"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Mortality_Rate_Difference, Unexposed_Rate, Exposed_Rate, Percent_Increase),
decimals = 2
)
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Rate Summary - Parasite Exposure | ||||
| Percent mortality and survival rates: A-T-P- vs A-T-P+ | ||||
| dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|
| A- T- P- | ||||
| 9 | 81 | 10.0 | 90 | 90.0 |
| A- T- P+ | ||||
| 8 | 82 | 8.9 | 90 | 91.1 |
| Mortality Rate Difference - Parasite Exposure | ||||
| Comparison of mortality rates between parasite-exposed and unexposed groups | ||||
| Comparison | Mortality_Rate_Difference | Unexposed_Rate | Exposed_Rate | Percent_Increase |
|---|---|---|---|---|
| A-T-P+ vs A-T-P- | −1.10 | 10.00 | 8.90 | −11.00 |
| Fisher's Exact Test Results - Parasite Exposure Mortality | |
| Comparison of A-T-P- vs A-T-P+ mortality rates | |
| Statistic | Value |
|---|---|
| Fisher's exact test p-value | 1.0000 |
| Odds ratio | 1.1381 |
| 95% CI lower | 0.3692 |
| 95% CI upper | 3.5739 |
| Chi-Square Test Results - Parasite Exposure Mortality | |
| Alternative test for A-T-P- vs A-T-P+ mortality rates | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 0 |
| Degrees of freedom | 1 |
| P-value | 1 |
Click on tabs to display tables. Scroll to see additional rows.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# Prepare data for statistical analysis - historical contingency mortality
historical_mortality_data <- data.mortality %>%
dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
dplyr::group_by(Treatment) %>%
dplyr::count() %>%
dplyr::mutate(
total_per_group = 90,
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1)
)
# Chi-square test for independence between treatment and mortality
set.seed(42)
historical_chi_square_test <- chisq.test(
matrix(c(historical_mortality_data$dead, historical_mortality_data$n),
ncol = 2, byrow = FALSE)
)
# Create mortality rate summary table for display
historical_mortality_rate_table <- historical_mortality_data %>%
dplyr::select(Treatment, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Summary - Historical Contingency",
subtitle = "Percent mortality and survival rates by treatment (parasite-exposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Chi-square test results
historical_chi_square_results <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(historical_chi_square_test$statistic, 4),
historical_chi_square_test$parameter,
round(historical_chi_square_test$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - Historical Contingency Mortality",
subtitle = "Test for independence between treatment and mortality (parasite-exposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
historical_pairwise_results <- list()
historical_treatments <- historical_mortality_data$Treatment
for(i in 1:(length(historical_treatments)-1)) {
for(j in (i+1):length(historical_treatments)) {
trt1 <- historical_treatments[i]
trt2 <- historical_treatments[j]
# Get data for these two treatments
data1 <- historical_mortality_data %>% dplyr::filter(Treatment == trt1)
data2 <- historical_mortality_data %>% dplyr::filter(Treatment == trt2)
# Create 2x2 contingency table
cont_table <- matrix(
c(data1$dead, data1$n, data2$dead, data2$n),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(trt1, trt2)
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
historical_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
Comparison = paste(trt1, "vs", trt2),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Combine pairwise results
historical_pairwise_summary <- do.call(rbind, historical_pairwise_results) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - Historical Contingency Mortality",
subtitle = "Post-hoc comparisons between treatment pairs (parasite-exposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
data.mortality %>%
dplyr::filter(Treatment %in% c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+")) %>%
dplyr::group_by(Treatment) %>%
dplyr::count() %>%
dplyr::mutate(
total_per_group = 90,
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1)
)
## # A tibble: 4 × 5
## # Groups: Treatment [4]
## Treatment n total_per_group dead percent_mortality
## <fct> <int> <dbl> <dbl> <dbl>
## 1 A- T- P+ 82 90 8 8.9
## 2 A+ T- P+ 80 90 10 11.1
## 3 A- T+ P+ 77 90 13 14.4
## 4 A+ T+ P+ 70 90 20 22.2
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Rate Summary - Historical Contingency | ||||
| Percent mortality and survival rates by treatment (parasite-exposed only) | ||||
| dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|
| A- T- P+ | ||||
| 8 | 82 | 8.9 | 90 | 91.1 |
| A+ T- P+ | ||||
| 10 | 80 | 11.1 | 90 | 88.9 |
| A- T+ P+ | ||||
| 13 | 77 | 14.4 | 90 | 85.6 |
| A+ T+ P+ | ||||
| 20 | 70 | 22.2 | 90 | 77.8 |
| Chi-Square Test Results - Historical Contingency Mortality | |
| Test for independence between treatment and mortality (parasite-exposed only) | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 7.561400 |
| Degrees of freedom | 3.000000 |
| P-value | 0.056002 |
| Pairwise Chi-Square Tests - Historical Contingency Mortality | |||
| Post-hoc comparisons between treatment pairs (parasite-exposed only) | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| A- T- P+ vs A+ T- P+ | 0.0617 | 0.803785 | No |
| A- T- P+ vs A- T+ P+ | 0.8625 | 0.353031 | No |
| A- T- P+ vs A+ T+ P+ | 5.1175 | 0.023686 | Yes |
| A+ T- P+ vs A- T+ P+ | 0.1994 | 0.655213 | No |
| A+ T- P+ vs A+ T+ P+ | 3.2400 | 0.071861 | No |
| A- T+ P+ vs A+ T+ P+ | 1.3358 | 0.247775 | No |
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# Statistical analysis for infection prevalence and worm burden
set.seed(42)
# 1. INFECTION PREVALENCE ANALYSIS
# Chi-square test for infection prevalence across treatments
infection_chi_square <- chisq.test(
data.infection %>%
dplyr::select(n_infected, n_fish) %>%
dplyr::mutate(n_uninfected = n_fish - n_infected) %>%
dplyr::select(n_infected, n_uninfected) %>%
as.matrix()
)
# Create infection prevalence contingency table
infection_contingency_table <- data.infection %>%
dplyr::mutate(
n_uninfected = n_fish - n_infected,
percent_uninfected = round((n_uninfected / n_fish) * 100, 1)
) %>%
dplyr::select(Treatment, n_infected, n_uninfected, n_fish, percent_infected, percent_uninfected) %>%
gt::gt() %>%
gt::tab_header(
title = "Infection Prevalence Contingency Table",
subtitle = "Number of infected vs uninfected fish by treatment"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(n_infected, n_uninfected, n_fish),
decimals = 0
)
# Chi-square test results for infection prevalence
infection_chi_results <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(infection_chi_square$statistic, 4),
infection_chi_square$parameter,
round(infection_chi_square$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - Infection Prevalence",
subtitle = "Test for independence between treatment and infection status"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise chi-square tests for infection prevalence
set.seed(42)
infection_pairwise_results <- list()
infection_treatments <- data.infection$Treatment
for(i in 1:(length(infection_treatments)-1)) {
for(j in (i+1):length(infection_treatments)) {
trt1 <- infection_treatments[i]
trt2 <- infection_treatments[j]
# Get data for these two treatments
data1 <- data.infection %>% dplyr::filter(Treatment == trt1)
data2 <- data.infection %>% dplyr::filter(Treatment == trt2)
# Create 2x2 contingency table
cont_table <- matrix(
c(data1$n_infected, data1$n_fish - data1$n_infected,
data2$n_infected, data2$n_fish - data2$n_infected),
nrow = 2, ncol = 2,
dimnames = list(
c("Infected", "Uninfected"),
c(trt1, trt2)
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
infection_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
Comparison = paste(trt1, "vs", trt2),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Combine pairwise results for infection prevalence
infection_pairwise_summary <- do.call(rbind, infection_pairwise_results) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - Infection Prevalence",
subtitle = "Post-hoc comparisons between treatment pairs"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# 2. WORM BURDEN ANALYSIS (TOTAL WORMS)
# Prepare individual-level data for negative binomial regression
set.seed(42)
individual_worm_data <- data.mortality %>%
dplyr::filter(Pathogen == 1) %>%
dplyr::select(Treatment, Total.Worm.Count) %>%
dplyr::filter(!is.na(Total.Worm.Count)) %>%
dplyr::mutate(
Treatment = factor(Treatment, levels = c("A- T- P+", "A+ T- P+", "A- T+ P+", "A+ T+ P+"))
)
# Fit negative binomial model for total worm counts
# Using MASS::glm.nb for negative binomial regression
library(MASS)
nb_model <- MASS::glm.nb(Total.Worm.Count ~ Treatment, data = individual_worm_data)
# Get model summary
nb_summary <- summary(nb_model)
nb_model.anova <- Anova(nb_model, type = 2)
# Extract coefficients and significance
nb_coefficients <- data.frame(
Treatment = c("Intercept", levels(individual_worm_data$Treatment)[-1]),
Estimate = round(nb_summary$coefficients[, 1], 4),
Std_Error = round(nb_summary$coefficients[, 2], 4),
z_value = round(nb_summary$coefficients[, 3], 4),
p_value = round(nb_summary$coefficients[, 4], 6)
) %>%
dplyr::mutate(
p_signif = ifelse(p_value < 0.001, "***",
ifelse(p_value < 0.01, "**",
ifelse(p_value < 0.05, "*", "ns")))
)
# Create worm burden summary table
worm_burden_data <- data.mortality %>%
dplyr::filter(Pathogen == 1) %>%
dplyr::group_by(Treatment) %>%
dplyr::summarise(
total_worms = sum(Total.Worm.Count, na.rm = TRUE),
mean_worms = mean(Total.Worm.Count, na.rm = TRUE),
median_worms = median(Total.Worm.Count, na.rm = TRUE),
sd_worms = sd(Total.Worm.Count, na.rm = TRUE),
n_fish = n(),
.groups = "drop"
) %>%
dplyr::mutate(
se_worms = sd_worms / sqrt(n_fish)
)
worm_burden_summary <- worm_burden_data %>%
dplyr::select(Treatment, total_worms, mean_worms, median_worms, sd_worms, se_worms, n_fish) %>%
gt::gt() %>%
gt::tab_header(
title = "Worm Burden Summary by Treatment",
subtitle = "Total, mean, median, standard deviation, and standard error of worm counts per treatment"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(total_worms, n_fish),
decimals = 0
) %>%
gt::fmt_number(
columns = c(mean_worms, median_worms, sd_worms, se_worms),
decimals = 2
)
# Negative binomial model results
nb_model_results <- nb_coefficients %>%
gt::gt() %>%
gt::tab_header(
title = "Negative Binomial Regression Results - Worm Burden",
subtitle = "Model: glm.nb(worm counts ~ Treatment)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Estimate, Std_Error, z_value, p_value),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
nb_model.anova %>%
broom::tidy() %>%
dplyr::mutate(
p_signif = ifelse(p.value < 0.001, "***",
ifelse(p.value < 0.01, "**",
ifelse(p.value < 0.05, "*", "ns")))
) %>%
gt::gt() %>%
gt::tab_header(
title = "Negative Binomial Regression ANOVA Results - Worm Burden",
subtitle = "Model: Anova(glm.nb(worm counts ~ Treatment))"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(p.value),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
| Negative Binomial Regression ANOVA Results - Worm Burden | ||||
| Model: Anova(glm.nb(worm counts ~ Treatment)) | ||||
| term | statistic | df | p.value | p_signif |
|---|---|---|---|---|
| Treatment | 9.044622 | 3 | 0.029 | * |
# Model fit statistics
nb_fit_stats <- data.frame(
Statistic = c("Log-likelihood", "AIC", "Theta", "SE Theta", "Deviance", "Residual Deviance"),
Value = c(
round(logLik(nb_model), 2),
round(AIC(nb_model), 2),
round(nb_model$theta, 4),
round(nb_model$SE.theta, 4),
round(nb_summary$deviance, 2),
round(sum(residuals(nb_model, type = "deviance")^2), 2)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Negative Binomial Model Fit Statistics",
subtitle = "Model diagnostics and fit measures"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise comparisons for worm burden using emmeans
set.seed(42)
library(emmeans)
# Get estimated marginal means and pairwise comparisons
nb_emmeans <- emmeans::emmeans(nb_model, ~ Treatment)
nb_pairwise <- pairs(nb_emmeans, adjust = "bonferroni")
# Convert to data frame for display
nb_pairwise_df <- as.data.frame(nb_pairwise) %>%
dplyr::mutate(
contrast = as.character(contrast),
estimate = round(estimate, 4),
SE = round(SE, 4),
z.ratio = round(z.ratio, 4),
p.value = round(p.value, 6),
p_signif = ifelse(p.value < 0.001, "***",
ifelse(p.value < 0.01, "**",
ifelse(p.value < 0.05, "*", "ns")))
) %>%
dplyr::select(contrast, estimate, SE, z.ratio, p.value, p_signif)
# Display pairwise comparisons for worm burden
nb_pairwise_summary <- nb_pairwise_df %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Comparisons - Worm Burden (Negative Binomial)",
subtitle = "Bonferroni-adjusted pairwise comparisons between treatments"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(estimate, SE, z.ratio, p.value),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
# Model diagnostics for negative binomial
# Check dispersion parameter and model adequacy
nb_diagnostics <- data.frame(
Statistic = c("Theta (dispersion parameter)", "SE Theta", "Deviance", "Degrees of Freedom", "Deviance/DF", "Model adequacy"),
Value = c(
round(nb_model$theta, 4),
round(nb_model$SE.theta, 4),
round(nb_summary$deviance, 2),
nb_summary$df.residual,
round(nb_summary$deviance / nb_summary$df.residual, 2),
ifelse(nb_summary$deviance / nb_summary$df.residual < 2, "Good fit", "Consider model diagnostics")
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Negative Binomial Model Diagnostics",
subtitle = "Checking model adequacy and dispersion"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Compare with Poisson model for reference
set.seed(42)
poisson_model <- glm(Total.Worm.Count ~ Treatment, data = individual_worm_data, family = "poisson")
poisson_summary <- summary(poisson_model)
# Model comparison
model_comparison <- data.frame(
Model = c("Poisson", "Negative Binomial"),
Log_Likelihood = c(round(logLik(poisson_model), 2), round(logLik(nb_model), 2)),
AIC = c(round(AIC(poisson_model), 2), round(AIC(nb_model), 2)),
Deviance = c(round(poisson_summary$deviance, 2), round(nb_summary$deviance, 2)),
DF = c(poisson_summary$df.residual, nb_summary$df.residual),
Deviance_DF_Ratio = c(
round(poisson_summary$deviance / poisson_summary$df.residual, 2),
round(nb_summary$deviance / nb_summary$df.residual, 2)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Model Comparison: Poisson vs Negative Binomial",
subtitle = "Comparing model fit and adequacy"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Log_Likelihood, AIC, Deviance, Deviance_DF_Ratio),
decimals = 2
) %>%
gt::fmt_number(
columns = DF,
decimals = 0
)
# 3. NORMALIZED WORM BURDEN ANALYSIS (WORMS PER INFECTED FISH)
# Calculate worms per infected fish to normalize for infection prevalence
# NOTE: This is ratio data, not count data, so we need different statistical approaches
set.seed(42)
# Create normalized worm burden summary table
normalized_worm_summary <- normalized_worm_data %>%
dplyr::select(Treatment, total_worms, n_infected, worms_per_infected, n_fish) %>%
dplyr::mutate(
infection_rate = round((n_infected / n_fish) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Normalized Worm Burden Summary by Treatment",
subtitle = "Worms per infected fish (normalized for infection prevalence)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(total_worms, n_infected, n_fish),
decimals = 0
) %>%
gt::fmt_number(
columns = c(worms_per_infected, infection_rate),
decimals = 2
)
# Create tank_normalized_for_test dataset for statistical analysis
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)
tank_normalized_for_test <- tank_normalized_data %>%
dplyr::filter(n_infected > 0) %>%
dplyr::mutate(
Treatment = factor(Treatment, levels = treatment_order)
)
# APPROPRIATE STATISTICAL ANALYSIS FOR RATIO DATA
# Method 1: One-way ANOVA on tank-level ratios (if assumptions are met)
set.seed(42)
# Use tank-level data for statistical testing
normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
normalized_anova_summary <- summary(normalized_anova)
# Extract ANOVA results using broom package for cleaner extraction
library(broom)
normalized_anova_tidy <- tidy(normalized_anova)
# Calculate F-statistic and p-value manually from ANOVA summary
f_stat <- normalized_anova_summary[[1]]$"Mean Sq"[1] / normalized_anova_summary[[1]]$"Mean Sq"[2]
p_val <- pf(f_stat, normalized_anova_summary[[1]]$"Df"[1], normalized_anova_summary[[1]]$"Df"[2], lower.tail = FALSE)
normalized_anova_results <- data.frame(
Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
Value = c(
round(f_stat, 4),
normalized_anova_summary[[1]]$"Df"[1],
normalized_anova_summary[[1]]$"Df"[2],
round(p_val, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "ANOVA Results - Normalized Worm Burden",
subtitle = "One-way ANOVA: worms per infected fish ~ Treatment (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc Tukey test for normalized worm burden
set.seed(42)
normalized_tukey <- TukeyHSD(normalized_anova)
# Convert Tukey results to data frame
normalized_tukey_df <- as.data.frame(normalized_tukey$Treatment) %>%
dplyr::mutate(
comparison = rownames(.),
diff = round(diff, 4),
lwr = round(lwr, 4),
upr = round(upr, 4),
p_adj = round(`p adj`, 6),
p_signif = ifelse(p_adj < 0.001, "***",
ifelse(p_adj < 0.01, "**",
ifelse(p_adj < 0.05, "*", "ns")))
) %>%
dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)
# Display Tukey test results
normalized_tukey_summary <- normalized_tukey_df %>%
gt::gt() %>%
gt::tab_header(
title = "Tukey HSD Test Results - Normalized Worm Burden",
subtitle = "Post-hoc pairwise comparisons for worms per infected fish (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(diff, lwr, upr, p_adj),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
# Method 2: Kruskal-Wallis test (non-parametric alternative)
set.seed(42)
normalized_kruskal <- kruskal.test(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
# Kruskal-Wallis results
normalized_kruskal_results <- data.frame(
Statistic = c("Kruskal-Wallis chi-squared", "Degrees of freedom", "P-value"),
Value = c(
round(normalized_kruskal$statistic, 4),
normalized_kruskal$parameter,
round(normalized_kruskal$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Kruskal-Wallis Test Results - Normalized Worm Burden",
subtitle = "Non-parametric test for differences in worms per infected fish across treatments (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Method 3: Permutation test (most appropriate for ratio data with small sample sizes)
set.seed(42)
library(coin)
# Create permutation test for normalized worm burden using tank-level data
normalized_permutation <- coin::oneway_test(
worms_per_infected ~ Treatment,
data = tank_normalized_for_test,
distribution = approximate(nresample = 10000)
)
# Permutation test results
normalized_permutation_results <- data.frame(
Statistic = c("Z statistic", "P-value (approximate)", "Number of permutations"),
Value = c(
round(statistic(normalized_permutation), 4),
round(pvalue(normalized_permutation), 6),
attr(pvalue(normalized_permutation), "nresample")
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Permutation Test Results - Normalized Worm Burden",
subtitle = "Non-parametric permutation test for ratio data (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Method 4: Bootstrap confidence intervals for ratios
set.seed(42)
library(boot)
##
## Attaching package: 'boot'
## The following object is masked from 'package:survival':
##
## aml
## The following object is masked from 'package:car':
##
## logit
# Function to calculate mean ratio for bootstrap
ratio_mean <- function(data, indices) {
d <- data[indices, ]
return(mean(d$worms_per_infected, na.rm = TRUE))
}
# Bootstrap analysis for each treatment using tank-level data
# We need multiple observations per treatment for bootstrap to work
bootstrap_results <- list()
for(trt in unique(tank_normalized_data$Treatment)) {
trt_data <- tank_normalized_data %>% dplyr::filter(Treatment == trt, n_infected > 0)
# Only proceed if we have multiple observations
if(nrow(trt_data) > 1) {
# Bootstrap with 10000 resamples
boot_result <- boot(trt_data, ratio_mean, R = 10000)
# Calculate confidence intervals
ci_result <- boot.ci(boot_result, type = "perc", conf = 0.95)
bootstrap_results[[trt]] <- data.frame(
Treatment = trt,
Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
Bootstrap_Mean = round(boot_result$t0, 3),
CI_Lower = round(ci_result$percent[4], 3),
CI_Upper = round(ci_result$percent[5], 3),
N_Tanks = nrow(trt_data)
)
} else {
# For treatments with only one observation, just report the mean
bootstrap_results[[trt]] <- data.frame(
Treatment = trt,
Mean_Ratio = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
Bootstrap_Mean = round(mean(trt_data$worms_per_infected, na.rm = TRUE), 3),
CI_Lower = NA,
CI_Upper = NA,
N_Tanks = nrow(trt_data)
)
}
}
# Combine bootstrap results
bootstrap_summary <- do.call(rbind, bootstrap_results) %>%
gt::gt() %>%
gt::tab_header(
title = "Bootstrap Confidence Intervals - Normalized Worm Burden",
subtitle = "95% confidence intervals for worms per infected fish by treatment (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Mean_Ratio, Bootstrap_Mean, CI_Lower, CI_Upper),
decimals = 3
)
# Method 5: Pairwise permutation tests
set.seed(42)
pairwise_permutation_results <- list()
treatments <- unique(tank_normalized_for_test$Treatment)
for(i in 1:(length(treatments)-1)) {
for(j in (i+1):length(treatments)) {
trt1 <- treatments[i]
trt2 <- treatments[j]
# Get data for these two treatments
data1 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt1)
data2 <- tank_normalized_for_test %>% dplyr::filter(Treatment == trt2)
# Combine data for permutation test
combined_data <- rbind(data1, data2)
# Perform permutation test
perm_test <- coin::oneway_test(
worms_per_infected ~ Treatment,
data = combined_data,
distribution = approximate(nresample = 10000)
)
pairwise_permutation_results[[paste(trt1, "vs", trt2)]] <- data.frame(
Comparison = paste(trt1, "vs", trt2),
Z_statistic = round(statistic(perm_test), 4),
P_value = round(pvalue(perm_test), 6),
Significant = ifelse(pvalue(perm_test) < 0.05, "Yes", "No")
)
}
}
# Combine pairwise permutation results
pairwise_permutation_summary <- do.call(rbind, pairwise_permutation_results) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Permutation Tests - Normalized Worm Burden",
subtitle = "Non-parametric pairwise comparisons for ratio data (tank-level data)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Model comparison and diagnostics
normalized_model_comparison <- data.frame(
Test = c("ANOVA", "Kruskal-Wallis", "Permutation Test"),
Statistic = c(
paste("F =", round(f_stat, 4)),
paste("χ² =", round(normalized_kruskal$statistic, 4)),
paste("Z =", round(statistic(normalized_permutation), 4))
),
P_value = c(
round(p_val, 6),
round(normalized_kruskal$p.value, 6),
round(pvalue(normalized_permutation), 6)
),
Assumptions = c(
"Normality, homogeneity of variance",
"None (non-parametric)",
"None (non-parametric)"
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Model Comparison - Normalized Worm Burden Analysis",
subtitle = "Comparing different statistical approaches for ratio data"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Tank-level normalized worm burden analysis
set.seed(42)
# Tank-level normalized worm burden summary
tank_normalized_summary <- tank_normalized_data %>%
dplyr::select(Treatment, Tank.ID, total_worms, n_infected, worms_per_infected, n_fish, infection_rate) %>%
gt::gt() %>%
gt::tab_header(
title = "Tank-Level Normalized Worm Burden Summary",
subtitle = "Worms per infected fish by treatment and tank"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(total_worms, n_infected, n_fish),
decimals = 0
) %>%
gt::fmt_number(
columns = c(worms_per_infected, infection_rate),
decimals = 2
)
# Statistical test for tank-level normalized worm burden
# Using ANOVA on worms per infected fish (excluding tanks with 0 infected fish)
set.seed(42)
# One-way ANOVA for worms per infected fish across treatments
tank_normalized_anova <- aov(worms_per_infected ~ Treatment, data = tank_normalized_for_test)
tank_normalized_anova_summary <- summary(tank_normalized_anova)
# Extract ANOVA results
tank_normalized_anova_results <- data.frame(
Statistic = c("F-value", "Degrees of freedom (Treatment)", "Degrees of freedom (Residuals)", "P-value"),
Value = c(
round(tank_normalized_anova_summary[[1]]$"F value"[1], 4),
tank_normalized_anova_summary[[1]]$"Df"[1],
tank_normalized_anova_summary[[1]]$"Df"[2],
round(tank_normalized_anova_summary[[1]]$"Pr(>F)"[1], 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "ANOVA Results - Tank-Level Normalized Worm Burden",
subtitle = "One-way ANOVA: worms per infected fish ~ Treatment"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc Tukey test for tank-level normalized worm burden
set.seed(42)
tank_normalized_tukey <- TukeyHSD(tank_normalized_anova)
# Convert Tukey results to data frame
tank_normalized_tukey_df <- as.data.frame(tank_normalized_tukey$Treatment) %>%
dplyr::mutate(
comparison = rownames(.),
diff = round(diff, 4),
lwr = round(lwr, 4),
upr = round(upr, 4),
p_adj = round(`p adj`, 6),
p_signif = ifelse(p_adj < 0.001, "***",
ifelse(p_adj < 0.01, "**",
ifelse(p_adj < 0.05, "*", "ns")))
) %>%
dplyr::select(comparison, diff, lwr, upr, p_adj, p_signif)
# Display Tukey test results
tank_normalized_tukey_summary <- tank_normalized_tukey_df %>%
gt::gt() %>%
gt::tab_header(
title = "Tukey HSD Test Results - Tank-Level Normalized Worm Burden",
subtitle = "Post-hoc pairwise comparisons for worms per infected fish"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(diff, lwr, upr, p_adj),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
Click on tabs to display tables. Scroll to see additional rows.
| Infection Prevalence Contingency Table | |||||
| Number of infected vs uninfected fish by treatment | |||||
| Treatment | n_infected | n_uninfected | n_fish | percent_infected | percent_uninfected |
|---|---|---|---|---|---|
| A- T- P+ | 31 | 51 | 82 | 37.8 | 62.2 |
| A+ T- P+ | 19 | 61 | 80 | 23.8 | 76.2 |
| A- T+ P+ | 21 | 56 | 77 | 27.3 | 72.7 |
| A+ T+ P+ | 15 | 55 | 70 | 21.4 | 78.6 |
| Chi-Square Test Results - Infection Prevalence | |
| Test for independence between treatment and infection status | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 6.16510 |
| Degrees of freedom | 3.00000 |
| P-value | 0.10385 |
| Pairwise Chi-Square Tests - Infection Prevalence | |||
| Post-hoc comparisons between treatment pairs | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| A- T- P+ vs A+ T- P+ | 3.1190 | 0.077384 | No |
| A- T- P+ vs A- T+ P+ | 1.5515 | 0.212910 | No |
| A- T- P+ vs A+ T+ P+ | 4.0541 | 0.044064 | Yes |
| A+ T- P+ vs A- T+ P+ | 0.1045 | 0.746535 | No |
| A+ T- P+ vs A+ T+ P+ | 0.0205 | 0.886027 | No |
| A- T+ P+ vs A+ T+ P+ | 0.3980 | 0.528098 | No |
| Worm Burden Summary by Treatment | ||||||
| Total, mean, median, standard deviation, and standard error of worm counts per treatment | ||||||
| Treatment | total_worms | mean_worms | median_worms | sd_worms | se_worms | n_fish |
|---|---|---|---|---|---|---|
| A- T- P+ | 167 | 2.04 | 0.00 | 5.16 | 0.57 | 82 |
| A+ T- P+ | 45 | 0.56 | 0.00 | 1.58 | 0.18 | 80 |
| A- T+ P+ | 121 | 1.57 | 0.00 | 3.59 | 0.41 | 77 |
| A+ T+ P+ | 77 | 1.10 | 0.00 | 2.82 | 0.34 | 70 |
| Negative Binomial Regression Results - Worm Burden | |||||
| Model: glm.nb(worm counts ~ Treatment) | |||||
| Treatment | Estimate | Std_Error | z_value | p_value | p_signif |
|---|---|---|---|---|---|
| Intercept | 0.711 | 0.292 | 2.433 | 0.015 | * |
| A+ T- P+ | −1.287 | 0.435 | −2.958 | 0.003 | ** |
| A- T+ P+ | −0.259 | 0.422 | −0.614 | 0.539 | ns |
| A+ T+ P+ | −0.616 | 0.438 | −1.407 | 0.159 | ns |
| Negative Binomial Model Fit Statistics | |
| Model diagnostics and fit measures | |
| Statistic | Value |
|---|---|
| Log-likelihood | -388.0400 |
| AIC | 786.0800 |
| Theta | 0.1534 |
| SE Theta | 0.0226 |
| Deviance | 189.0100 |
| Residual Deviance | 189.0100 |
| Pairwise Comparisons - Worm Burden (Negative Binomial) | |||||
| Bonferroni-adjusted pairwise comparisons between treatments | |||||
| contrast | estimate | SE | z.ratio | p.value | p_signif |
|---|---|---|---|---|---|
| (A- T- P+) - (A+ T- P+) | 1.287 | 0.435 | 2.958 | 0.019 | * |
| (A- T- P+) - (A- T+ P+) | 0.259 | 0.422 | 0.614 | 1.000 | ns |
| (A- T- P+) - (A+ T+ P+) | 0.616 | 0.438 | 1.407 | 0.956 | ns |
| (A+ T- P+) - (A- T+ P+) | −1.027 | 0.443 | −2.317 | 0.123 | ns |
| (A+ T- P+) - (A+ T+ P+) | −0.671 | 0.458 | −1.464 | 0.859 | ns |
| (A- T+ P+) - (A+ T+ P+) | 0.357 | 0.446 | 0.800 | 1.000 | ns |
| Negative Binomial Model Diagnostics | |
| Checking model adequacy and dispersion | |
| Statistic | Value |
|---|---|
| Theta (dispersion parameter) | 0.1534 |
| SE Theta | 0.0226 |
| Deviance | 189.01 |
| Degrees of Freedom | 305 |
| Deviance/DF | 0.62 |
| Model adequacy | Good fit |
| Model Comparison: Poisson vs Negative Binomial | |||||
| Comparing model fit and adequacy | |||||
| Model | Log_Likelihood | AIC | Deviance | DF | Deviance_DF_Ratio |
|---|---|---|---|---|---|
| Poisson | −800.23 | 1,608.45 | 1,338.72 | 305 | 4.39 |
| Negative Binomial | −388.04 | 786.08 | 189.01 | 305 | 0.62 |
| Normalized Worm Burden Summary by Treatment | |||||
| Worms per infected fish (normalized for infection prevalence) | |||||
| Treatment | total_worms | n_infected | worms_per_infected | n_fish | infection_rate |
|---|---|---|---|---|---|
| A- T- P+ | 167 | 31 | 5.39 | 82 | 37.80 |
| A+ T- P+ | 45 | 19 | 2.37 | 80 | 23.80 |
| A- T+ P+ | 121 | 21 | 5.76 | 77 | 27.30 |
| A+ T+ P+ | 77 | 15 | 5.13 | 70 | 21.40 |
| ANOVA Results - Normalized Worm Burden | |
| One-way ANOVA: worms per infected fish ~ Treatment (tank-level data) | |
| Statistic | Value |
|---|---|
| F-value | 1.857900 |
| Degrees of freedom (Treatment) | 3.000000 |
| Degrees of freedom (Residuals) | 8.000000 |
| P-value | 0.215039 |
| Tukey HSD Test Results - Normalized Worm Burden | |||||
| Post-hoc pairwise comparisons for worms per infected fish (tank-level data) | |||||
| comparison | diff | lwr | upr | p_adj | p_signif |
|---|---|---|---|---|---|
| A+ T- P+-A- T- P+ | −3.012 | −8.471 | 2.448 | 0.354 | ns |
| A- T+ P+-A- T- P+ | 0.198 | −5.262 | 5.657 | 0.999 | ns |
| A+ T+ P+-A- T- P+ | 0.536 | −4.924 | 5.995 | 0.988 | ns |
| A- T+ P+-A+ T- P+ | 3.210 | −2.250 | 8.669 | 0.307 | ns |
| A+ T+ P+-A+ T- P+ | 3.548 | −1.912 | 9.007 | 0.238 | ns |
| A+ T+ P+-A- T+ P+ | 0.338 | −5.122 | 5.798 | 0.997 | ns |
| Kruskal-Wallis Test Results - Normalized Worm Burden | |
| Non-parametric test for differences in worms per infected fish across treatments (tank-level data) | |
| Statistic | Value |
|---|---|
| Kruskal-Wallis chi-squared | 3.769200 |
| Degrees of freedom | 3.000000 |
| P-value | 0.287485 |
| Permutation Test Results - Normalized Worm Burden | |
| Non-parametric permutation test for ratio data (tank-level data) | |
| Statistic | Value |
|---|---|
| Z statistic | 4.5168 |
| P-value (approximate) | 0.2356 |
| Number of permutations | 10000.0000 |
| Bootstrap Confidence Intervals - Normalized Worm Burden | |||||
| 95% confidence intervals for worms per infected fish by treatment (tank-level data) | |||||
| Treatment | Mean_Ratio | Bootstrap_Mean | CI_Lower | CI_Upper | N_Tanks |
|---|---|---|---|---|---|
| A- T- P+ | 5.095 | 5.095 | 3.286 | 6.167 | 3 |
| A+ T- P+ | 2.083 | 2.083 | 1.000 | 3.625 | 3 |
| A- T+ P+ | 5.293 | 5.293 | 2.333 | 7.545 | 3 |
| A+ T+ P+ | 5.631 | 5.631 | 3.143 | 8.000 | 3 |
| Pairwise Permutation Tests - Normalized Worm Burden | |||
| Non-parametric pairwise comparisons for ratio data (tank-level data) | |||
| Comparison | Z_statistic | P_value | Significant |
|---|---|---|---|
| A- T- P+ vs A+ T- P+ | 1.7453 | 0.2013 | No |
| A- T- P+ vs A- T+ P+ | -0.1230 | 0.9007 | No |
| A- T- P+ vs A+ T+ P+ | -0.3536 | 1.0000 | No |
| A+ T- P+ vs A- T+ P+ | -1.5176 | 0.2005 | No |
| A+ T- P+ vs A+ T+ P+ | -1.6551 | 0.2003 | No |
| A- T+ P+ vs A+ T+ P+ | -0.1804 | 0.8006 | No |
| Model Comparison - Normalized Worm Burden Analysis | |||
| Comparing different statistical approaches for ratio data | |||
| Test | Statistic | P_value | Assumptions |
|---|---|---|---|
| ANOVA | F = 1.8579 | 0.215039 | Normality, homogeneity of variance |
| Kruskal-Wallis | χ² = 3.7692 | 0.287485 | None (non-parametric) |
| Permutation Test | Z = 4.5168 | 0.235600 | None (non-parametric) |
| Tank-Level Normalized Worm Burden Summary | ||||||
| Worms per infected fish by treatment and tank | ||||||
| Treatment | Tank.ID | total_worms | n_infected | worms_per_infected | n_fish | infection_rate |
|---|---|---|---|---|---|---|
| A- T- P+ | 4 | 70 | 12 | 5.83 | 29 | 41.40 |
| A- T- P+ | 5 | 23 | 7 | 3.29 | 25 | 28.00 |
| A- T- P+ | 6 | 74 | 12 | 6.17 | 28 | 42.90 |
| A+ T- P+ | 16 | 13 | 8 | 1.62 | 29 | 27.60 |
| A+ T- P+ | 17 | 29 | 8 | 3.62 | 28 | 28.60 |
| A+ T- P+ | 18 | 3 | 3 | 1.00 | 23 | 13.00 |
| A- T+ P+ | 10 | 24 | 4 | 6.00 | 21 | 19.00 |
| A- T+ P+ | 11 | 14 | 6 | 2.33 | 30 | 20.00 |
| A- T+ P+ | 12 | 83 | 11 | 7.55 | 26 | 42.30 |
| A+ T+ P+ | 22 | 22 | 7 | 3.14 | 28 | 25.00 |
| A+ T+ P+ | 23 | 23 | 4 | 5.75 | 23 | 17.40 |
| A+ T+ P+ | 24 | 32 | 4 | 8.00 | 19 | 21.10 |
| ANOVA Results - Tank-Level Normalized Worm Burden | |
| One-way ANOVA: worms per infected fish ~ Treatment | |
| Statistic | Value |
|---|---|
| F-value | 1.857900 |
| Degrees of freedom (Treatment) | 3.000000 |
| Degrees of freedom (Residuals) | 8.000000 |
| P-value | 0.215039 |
| Tukey HSD Test Results - Tank-Level Normalized Worm Burden | |||||
| Post-hoc pairwise comparisons for worms per infected fish | |||||
| comparison | diff | lwr | upr | p_adj | p_signif |
|---|---|---|---|---|---|
| A+ T- P+-A- T- P+ | −3.012 | −8.471 | 2.448 | 0.354 | ns |
| A- T+ P+-A- T- P+ | 0.198 | −5.262 | 5.657 | 0.999 | ns |
| A+ T+ P+-A- T- P+ | 0.536 | −4.924 | 5.995 | 0.988 | ns |
| A- T+ P+-A+ T- P+ | 3.210 | −2.250 | 8.669 | 0.307 | ns |
| A+ T+ P+-A+ T- P+ | 3.548 | −1.912 | 9.007 | 0.238 | ns |
| A+ T+ P+-A- T+ P+ | 0.338 | −5.122 | 5.798 | 0.997 | ns |
Click on tabs to display tables. Scroll to see additional rows.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's colour values.
# Prepare data for statistical analysis - recovery mortality
recovery_mortality_data <- data.mortality %>%
dplyr::filter(Treatment %in% c("A- T- P-", "A+ T- P-", "A- T+ P-", "A+ T+ P-")) %>%
dplyr::group_by(Treatment) %>%
dplyr::count() %>%
dplyr::mutate(
total_per_group = 90,
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1)
)
# Chi-square test for independence between treatment and mortality
set.seed(42)
recovery_chi_square_test <- chisq.test(
matrix(c(recovery_mortality_data$dead, recovery_mortality_data$n),
ncol = 2, byrow = FALSE)
)
# Create mortality rate summary table for display
recovery_mortality_rate_table <- recovery_mortality_data %>%
dplyr::select(Treatment, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Summary - Recovery",
subtitle = "Percent mortality and survival rates by treatment (parasite-unexposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Chi-square test results
recovery_chi_square_results <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(recovery_chi_square_test$statistic, 4),
recovery_chi_square_test$parameter,
round(recovery_chi_square_test$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - Recovery Mortality",
subtitle = "Test for independence between treatment and mortality (parasite-unexposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
recovery_pairwise_results <- list()
recovery_treatments <- recovery_mortality_data$Treatment
for(i in 1:(length(recovery_treatments)-1)) {
for(j in (i+1):length(recovery_treatments)) {
trt1 <- recovery_treatments[i]
trt2 <- recovery_treatments[j]
# Get data for these two treatments
data1 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt1)
data2 <- recovery_mortality_data %>% dplyr::filter(Treatment == trt2)
# Create 2x2 contingency table
cont_table <- matrix(
c(data1$dead, data1$n, data2$dead, data2$n),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(trt1, trt2)
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
recovery_pairwise_results[[paste(trt1, "vs", trt2)]] <- data.frame(
Comparison = paste(trt1, "vs", trt2),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Combine pairwise results
recovery_pairwise_summary <- do.call(rbind, recovery_pairwise_results) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - Recovery Mortality",
subtitle = "Post-hoc comparisons between treatment pairs (parasite-unexposed only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Rate Summary - Recovery | ||||
| Percent mortality and survival rates by treatment (parasite-unexposed only) | ||||
| dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|
| A- T- P- | ||||
| 9 | 81 | 10.0 | 90 | 90.0 |
| A+ T- P- | ||||
| 14 | 76 | 15.6 | 90 | 84.4 |
| A- T+ P- | ||||
| 19 | 71 | 21.1 | 90 | 78.9 |
| A+ T+ P- | ||||
| 30 | 60 | 33.3 | 90 | 66.7 |
| Chi-Square Test Results - Recovery Mortality | |
| Test for independence between treatment and mortality (parasite-unexposed only) | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 16.805600 |
| Degrees of freedom | 3.000000 |
| P-value | 0.000775 |
| Pairwise Chi-Square Tests - Recovery Mortality | |||
| Post-hoc comparisons between treatment pairs (parasite-unexposed only) | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| A- T- P- vs A+ T- P- | 0.7976 | 0.371823 | No |
| A- T- P- vs A- T+ P- | 3.4258 | 0.064187 | No |
| A- T- P- vs A+ T+ P- | 13.0933 | 0.000296 | Yes |
| A+ T- P- vs A- T+ P- | 0.5937 | 0.440995 | No |
| A+ T- P- vs A+ T+ P- | 6.7680 | 0.009280 | Yes |
| A- T+ P- vs A+ T+ P- | 2.8042 | 0.094019 | No |
Click on tabs to display tables. Scroll to see additional rows.
# Prepare data for statistical analysis - percent mortality rates by History
history_mortality_stats_data <- data.mortality %>%
dplyr::group_by(History) %>%
dplyr::count() %>%
dplyr::mutate(
# Calculate expected totals: 90 fish per treatment group
total_per_group = dplyr::case_when(
History == 0 ~ 90 * 2, # 2 treatment groups: A- T- P- and A- T- P+
History == 1 ~ 90 * 4, # 4 treatment groups: A+ T- P-, A+ T- P+, A- T+ P-, A- T+ P+
History == 2 ~ 90 * 2, # 2 treatment groups: A+ T+ P- and A+ T+ P+
TRUE ~ 0
),
dead = total_per_group - n,
percent_mortality = round((dead / total_per_group) * 100, 1),
History_Label = dplyr::case_when(
History == 0 ~ "No prior stressors",
History == 1 ~ "One prior stressor",
History == 2 ~ "Two prior stressors",
TRUE ~ "Unknown"
)
)
# Chi-square test for independence between history and mortality
set.seed(42)
history_chi_square_test <- chisq.test(
matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n),
ncol = 2, byrow = FALSE)
)
# Create mortality rate summary table for display
history_mortality_rate_table <- history_mortality_stats_data %>%
dplyr::select(History, History_Label, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Rate Summary - History of Stressors",
subtitle = "Percent mortality and survival rates by history of stressors"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Chi-square test results
history_test_stat_table <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(history_chi_square_test$statistic, 4),
history_chi_square_test$parameter,
round(history_chi_square_test$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - History of Stressors Mortality",
subtitle = "Test for independence between history of stressors and mortality"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Post-hoc pairwise comparisons using chi-square tests
set.seed(42)
history_pairwise_results <- list()
history_levels <- history_mortality_stats_data$History
for(i in 1:(length(history_levels)-1)) {
for(j in (i+1):length(history_levels)) {
hist1 <- history_levels[i]
hist2 <- history_levels[j]
# Get data for these two history levels
data1 <- history_mortality_stats_data %>% dplyr::filter(History == hist1)
data2 <- history_mortality_stats_data %>% dplyr::filter(History == hist2)
# Create 2x2 contingency table
cont_table <- matrix(
c(data1$dead, data1$n, data2$dead, data2$n),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(paste("History", hist1), paste("History", hist2))
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
history_pairwise_results[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
Comparison = paste("History", hist1, "vs", "History", hist2),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Combine pairwise results
history_pairwise_summary <- do.call(rbind, history_pairwise_results) %>%
dplyr::arrange(desc(Significant), P_value) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - History of Stressors Mortality",
subtitle = "Post-hoc comparisons between history levels"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Show only significant pairwise comparisons
history_pairwise_summary__SigOnly <- do.call(rbind, history_pairwise_results) %>%
dplyr::arrange(desc(Significant), P_value) %>%
dplyr::filter(Significant == "Yes") %>%
gt::gt() %>%
gt::tab_header(
title = "Significant Pairwise Chi-Square Tests - History of Stressors Mortality",
subtitle = "Post-hoc comparisons between history levels (significant only)"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Additional analysis: Trend analysis across history levels
# Calculate trend in mortality rates
history_trend_data <- history_mortality_stats_data %>%
dplyr::arrange(History) %>%
dplyr::select(History, percent_mortality)
# Linear trend test (Cochran-Armitage trend test equivalent using chi-square)
set.seed(42)
# Create trend weights (0, 1, 2 for history levels)
trend_weights <- c(0, 1, 2)
trend_chi_square <- chisq.test(
matrix(c(history_mortality_stats_data$dead, history_mortality_stats_data$n),
ncol = 2, byrow = FALSE),
simulate.p.value = TRUE
)
# Calculate correlation between history level and mortality rate
history_correlation <- cor.test(
history_trend_data$History,
history_trend_data$percent_mortality,
method = "pearson"
)
# Trend analysis results
history_trend_results <- data.frame(
Statistic = c("Correlation coefficient (r)", "P-value (correlation)", "Chi-square trend test p-value"),
Value = c(
round(history_correlation$estimate, 4),
round(history_correlation$p.value, 6),
round(trend_chi_square$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Trend Analysis - History of Stressors Mortality",
subtitle = "Testing for linear trend in mortality rates across history levels"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Effect size analysis: Calculate odds ratios between history levels
history_effect_sizes <- data.frame(
Comparison = c("History 1 vs History 0", "History 2 vs History 0", "History 2 vs History 1"),
Odds_Ratio = c(
round((history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[1]) /
(history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[2]), 3),
round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[1]) /
(history_mortality_stats_data$dead[1] * history_mortality_stats_data$n[3]), 3),
round((history_mortality_stats_data$dead[3] * history_mortality_stats_data$n[2]) /
(history_mortality_stats_data$dead[2] * history_mortality_stats_data$n[3]), 3)
),
Mortality_Rate_1 = c(
round(history_mortality_stats_data$percent_mortality[2], 1),
round(history_mortality_stats_data$percent_mortality[3], 1),
round(history_mortality_stats_data$percent_mortality[3], 1)
),
Mortality_Rate_2 = c(
round(history_mortality_stats_data$percent_mortality[1], 1),
round(history_mortality_stats_data$percent_mortality[1], 1),
round(history_mortality_stats_data$percent_mortality[2], 1)
)
) %>%
dplyr::mutate(
Interpretation = dplyr::case_when(
Odds_Ratio > 1 ~ "Higher mortality in first group",
Odds_Ratio < 1 ~ "Lower mortality in first group",
Odds_Ratio == 1 ~ "No difference"
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Effect Size Analysis - History of Stressors Mortality",
subtitle = "Odds ratios comparing mortality rates between history levels"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Odds_Ratio, Mortality_Rate_1, Mortality_Rate_2),
decimals = 1
)
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Rate Summary - History of Stressors | |||||
| Percent mortality and survival rates by history of stressors | |||||
| History_Label | dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|---|
| 0 | |||||
| No prior stressors | 17 | 163 | 9.4 | 180 | 90.6 |
| 1 | |||||
| One prior stressor | 56 | 304 | 15.6 | 360 | 84.4 |
| 2 | |||||
| Two prior stressors | 50 | 130 | 27.8 | 180 | 72.2 |
| Chi-Square Test Results - History of Stressors Mortality | |
| Test for independence between history of stressors and mortality | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 22.542000 |
| Degrees of freedom | 2.000000 |
| P-value | 0.000013 |
| Pairwise Chi-Square Tests - History of Stressors Mortality | |||
| Post-hoc comparisons between history levels | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| History 0 vs History 2 | 18.7785 | 0.000015 | Yes |
| History 1 vs History 2 | 10.6010 | 0.001130 | Yes |
| History 0 vs History 1 | 3.3284 | 0.068094 | No |
| Significant Pairwise Chi-Square Tests - History of Stressors Mortality | |||
| Post-hoc comparisons between history levels (significant only) | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| History 0 vs History 2 | 18.7785 | 0.000015 | Yes |
| History 1 vs History 2 | 10.6010 | 0.001130 | Yes |
| Trend Analysis - History of Stressors Mortality | |
| Testing for linear trend in mortality rates across history levels | |
| Statistic | Value |
|---|---|
| Correlation coefficient (r) | 0.982700 |
| P-value (correlation) | 0.118467 |
| Chi-square trend test p-value | 0.000500 |
| Effect Size Analysis - History of Stressors Mortality | ||||
| Odds ratios comparing mortality rates between history levels | ||||
| Comparison | Odds_Ratio | Mortality_Rate_1 | Mortality_Rate_2 | Interpretation |
|---|---|---|---|---|
| History 1 vs History 0 | 1.8 | 15.6 | 9.4 | Higher mortality in first group |
| History 2 vs History 0 | 3.7 | 27.8 | 9.4 | Higher mortality in first group |
| History 2 vs History 1 | 2.1 | 27.8 | 15.6 | Higher mortality in first group |
Click on tabs to display tables. Scroll to see additional rows.
# Prepare data for statistical analysis - History × Pathogen interaction
interaction_stats_data <- history_pathogen_interaction_data %>%
dplyr::select(History, Pathogen, History_Label, Pathogen_Label, dead, n, percent_mortality) %>%
dplyr::rename(Survived = n) %>%
dplyr::mutate(
Total = dead + Survived,
percent_survived = round((Survived / Total) * 100, 1)
)
# 1. Chi-square test for independence between History and Pathogen exposure
set.seed(42)
# Create 3x2 contingency table (History × Pathogen)
interaction_contingency <- matrix(
c(interaction_stats_data$dead),
nrow = 3, ncol = 2,
dimnames = list(
c("History 0", "History 1", "History 2"),
c("Unexposed", "Exposed")
)
)
interaction_chi_square <- chisq.test(interaction_contingency)
# 2. Test for interaction effect using logistic regression approach
# Create individual-level data for interaction analysis
set.seed(42)
individual_interaction_data <- data.mortality %>%
dplyr::select(History, Pathogen, Sample) %>%
dplyr::mutate(
# Create binary outcome: 1 = survived, 0 = died (assuming all sampled fish survived)
Outcome = 1, # All fish in the dataset survived (were sampled)
History = factor(History, levels = c(0, 1, 2)),
Pathogen = factor(Pathogen, levels = c(0, 1))
)
# Fit logistic regression model with interaction
interaction_model <- glm(Outcome ~ History * Pathogen,
data = individual_interaction_data,
family = "binomial")
## Warning: glm.fit: algorithm did not converge
# Get model summary
interaction_summary <- summary(interaction_model)
# Extract coefficients and significance
interaction_coefficients <- data.frame(
Term = rownames(interaction_summary$coefficients),
Estimate = round(interaction_summary$coefficients[, 1], 4),
Std_Error = round(interaction_summary$coefficients[, 2], 4),
z_value = round(interaction_summary$coefficients[, 3], 4),
p_value = round(interaction_summary$coefficients[, 4], 6)
) %>%
dplyr::mutate(
p_signif = ifelse(p_value < 0.001, "***",
ifelse(p_value < 0.01, "**",
ifelse(p_value < 0.05, "*", "ns")))
)
# 3. Create summary tables for display
interaction_summary_table <- interaction_stats_data %>%
gt::gt() %>%
gt::tab_header(
title = "Mortality Summary - History × Pathogen Interaction",
subtitle = "Percent mortality and survival rates by history and pathogen exposure"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(dead, Survived, Total),
decimals = 0
) %>%
gt::fmt_number(
columns = c(percent_mortality, percent_survived),
decimals = 1
)
# Chi-square test results
interaction_chi_results <- data.frame(
Statistic = c("Chi-square statistic", "Degrees of freedom", "P-value"),
Value = c(
round(interaction_chi_square$statistic, 4),
interaction_chi_square$parameter,
round(interaction_chi_square$p.value, 6)
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Chi-Square Test Results - History × Pathogen Interaction",
subtitle = "Test for independence between history and pathogen exposure"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# Logistic regression results
interaction_model_results <- interaction_coefficients %>%
gt::gt() %>%
gt::tab_header(
title = "Logistic Regression Results - History × Pathogen Interaction",
subtitle = "Model: Outcome ~ History * Pathogen"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Estimate, Std_Error, z_value, p_value),
decimals = 3
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = p_signif,
rows = p_signif != "ns"
)
)
# 4. Calculate effect sizes for each comparison
# Calculate odds ratios for each history level comparing exposed vs unexposed
set.seed(42)
effect_size_comparisons <- data.frame(
Comparison = c(
"History 0: Exposed vs Unexposed",
"History 1: Exposed vs Unexposed",
"History 2: Exposed vs Unexposed"
),
Odds_Ratio = c(
# History 0: Exposed vs Unexposed
round((interaction_stats_data$dead[2] * interaction_stats_data$Survived[1]) /
(interaction_stats_data$dead[1] * interaction_stats_data$Survived[2]), 3),
# History 1: Exposed vs Unexposed
round((interaction_stats_data$dead[4] * interaction_stats_data$Survived[3]) /
(interaction_stats_data$dead[3] * interaction_stats_data$Survived[4]), 3),
# History 2: Exposed vs Unexposed
round((interaction_stats_data$dead[6] * interaction_stats_data$Survived[5]) /
(interaction_stats_data$dead[5] * interaction_stats_data$Survived[6]), 3)
),
Mortality_Exposed = c(
round(interaction_stats_data$percent_mortality[2], 1),
round(interaction_stats_data$percent_mortality[4], 1),
round(interaction_stats_data$percent_mortality[6], 1)
),
Mortality_Unexposed = c(
round(interaction_stats_data$percent_mortality[1], 1),
round(interaction_stats_data$percent_mortality[3], 1),
round(interaction_stats_data$percent_mortality[5], 1)
)
) %>%
dplyr::mutate(
Mortality_Difference = Mortality_Exposed - Mortality_Unexposed,
Interpretation = dplyr::case_when(
Odds_Ratio > 1 ~ "Higher mortality when exposed",
Odds_Ratio < 1 ~ "Lower mortality when exposed",
Odds_Ratio == 1 ~ "No difference"
)
)
effect_size_comparisons_table <- effect_size_comparisons %>%
gt::gt() %>%
gt::tab_header(
title = "Effect Size Analysis - History × Pathogen Interaction",
subtitle = "Odds ratios comparing exposed vs unexposed within each history level"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Odds_Ratio, Mortality_Exposed, Mortality_Unexposed, Mortality_Difference),
decimals = 1
)
# 5. Test for trend in pathogen effect across history levels
set.seed(42)
# Calculate correlation between history level and the difference in mortality (exposed - unexposed)
trend_data <- effect_size_comparisons %>%
dplyr::mutate(History_Level = c(0, 1, 2))
pathogen_effect_trend <- cor.test(
trend_data$History_Level,
trend_data$Mortality_Difference,
method = "pearson"
)
# Trend analysis results
pathogen_effect_trend_results <- data.frame(
Statistic = c("Correlation coefficient (r)", "P-value", "Interpretation"),
Value = c(
round(pathogen_effect_trend$estimate, 4),
round(pathogen_effect_trend$p.value, 6),
ifelse(pathogen_effect_trend$p.value < 0.05,
"Significant trend in pathogen effect across history levels",
"No significant trend in pathogen effect across history levels")
)
) %>%
gt::gt() %>%
gt::tab_header(
title = "Trend Analysis - Pathogen Effect Across History Levels",
subtitle = "Testing whether the effect of pathogen exposure changes across history levels"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
)
# 6. Post-hoc pairwise comparisons within each pathogen group
set.seed(42)
# Within unexposed group
unexposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 0)
unexposed_pairwise <- list()
for(i in 1:(nrow(unexposed_data)-1)) {
for(j in (i+1):nrow(unexposed_data)) {
hist1 <- unexposed_data$History[i]
hist2 <- unexposed_data$History[j]
# Create 2x2 contingency table
cont_table <- matrix(
c(unexposed_data$dead[i], unexposed_data$Survived[i],
unexposed_data$dead[j], unexposed_data$Survived[j]),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(paste("History", hist1), paste("History", hist2))
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
unexposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
Comparison = paste("History", hist1, "vs", "History", hist2, "(Unexposed)"),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# Within exposed group
exposed_data <- interaction_stats_data %>% dplyr::filter(Pathogen == 1)
exposed_pairwise <- list()
for(i in 1:(nrow(exposed_data)-1)) {
for(j in (i+1):nrow(exposed_data)) {
hist1 <- exposed_data$History[i]
hist2 <- exposed_data$History[j]
# Create 2x2 contingency table
cont_table <- matrix(
c(exposed_data$dead[i], exposed_data$Survived[i],
exposed_data$dead[j], exposed_data$Survived[j]),
nrow = 2, ncol = 2,
dimnames = list(
c("Dead", "Survived"),
c(paste("History", hist1), paste("History", hist2))
)
)
# Perform chi-square test
test_result <- chisq.test(cont_table)
exposed_pairwise[[paste("History", hist1, "vs", "History", hist2)]] <- data.frame(
Comparison = paste("History", hist1, "vs", "History", hist2, "(Exposed)"),
Chi_square = round(test_result$statistic, 4),
P_value = round(test_result$p.value, 6),
Significant = ifelse(test_result$p.value < 0.05, "Yes", "No")
)
}
}
# 7. Pairwise comparisons between same history types across pathogen exposures
set.seed(42)
same_history_pathogen_pairwise <- list()
# For each history level, compare unexposed vs exposed
for(hist_level in c(0, 1, 2)) {
# Get data for this history level
unexposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 0)
exposed_row <- interaction_stats_data %>% dplyr::filter(History == hist_level & Pathogen == 1)
if(nrow(unexposed_row) > 0 && nrow(exposed_row) > 0) {
# Calculate percent mortality for each group
unexposed_mortality <- unexposed_row$percent_mortality
exposed_mortality <- exposed_row$percent_mortality
# Calculate the difference in percent mortality
mortality_difference <- exposed_mortality - unexposed_mortality
# Perform t-test on percent mortality rates
# Since we have single values, we'll use a simple comparison
# For small sample sizes, we can use a z-test for proportions
# Calculate standard error of the difference
unexposed_se <- sqrt((unexposed_mortality * (100 - unexposed_mortality)) / unexposed_row$Total)
exposed_se <- sqrt((exposed_mortality * (100 - exposed_mortality)) / exposed_row$Total)
se_diff <- sqrt(unexposed_se^2 + exposed_se^2)
# Calculate z-statistic
z_stat <- mortality_difference / se_diff
p_value <- 2 * (1 - pnorm(abs(z_stat)))
same_history_pathogen_pairwise[[paste("History", hist_level)]] <- data.frame(
Comparison = paste("History", hist_level, "(Unexposed) vs History", hist_level, "(Exposed)"),
Unexposed_Mortality = round(unexposed_mortality, 1),
Exposed_Mortality = round(exposed_mortality, 1),
Mortality_Difference = round(mortality_difference, 1),
Z_statistic = round(z_stat, 4),
P_value = round(p_value, 6),
Significant = ifelse(p_value < 0.05, "Yes", "No")
)
}
}
# Create separate tables for different types of pairwise comparisons
# Table 1: Within unexposed group comparisons
unexposed_pairwise_results <- do.call(rbind, unexposed_pairwise) %>%
dplyr::arrange(desc(Significant), P_value) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - Within Unexposed Group",
subtitle = "Comparing history levels within unexposed fish"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Table 2: Within exposed group comparisons
exposed_pairwise_results <- do.call(rbind, exposed_pairwise) %>%
dplyr::arrange(desc(Significant), P_value) %>%
gt::gt() %>%
gt::tab_header(
title = "Pairwise Chi-Square Tests - Within Exposed Group",
subtitle = "Comparing history levels within exposed fish"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
# Create separate table for same history comparisons across pathogen exposures
same_history_pathogen_results <- do.call(rbind, same_history_pathogen_pairwise) %>%
dplyr::arrange(desc(Significant), P_value) %>%
gt::gt() %>%
gt::tab_header(
title = "Pathogen Effect Within Each History Level",
subtitle = "Comparing percent mortality between unexposed vs exposed within the same history of stressors"
) %>%
gt::tab_style(
style = gt::cell_fill(color = "#f7f7f7"),
locations = gt::cells_column_labels()
) %>%
gt::fmt_number(
columns = c(Unexposed_Mortality, Exposed_Mortality, Mortality_Difference),
decimals = 1
) %>%
gt::fmt_number(
columns = c(Z_statistic, P_value),
decimals = 4
) %>%
gt::tab_style(
style = gt::cell_fill(color = "lightgreen"),
locations = gt::cells_body(
columns = Significant,
rows = Significant == "Yes"
)
)
Click on tabs to display tables. Scroll to see additional rows.
| Mortality Summary - History × Pathogen Interaction | ||||||
| Percent mortality and survival rates by history and pathogen exposure | ||||||
| History_Label | Pathogen_Label | dead | Survived | percent_mortality | Total | percent_survived |
|---|---|---|---|---|---|---|
| 0 - 0 | ||||||
| No prior stressors | Unexposed | 9 | 81 | 10.0 | 90 | 90.0 |
| 0 - 1 | ||||||
| No prior stressors | Exposed | 8 | 82 | 8.9 | 90 | 91.1 |
| 1 - 0 | ||||||
| One prior stressor | Unexposed | 33 | 147 | 18.3 | 180 | 81.7 |
| 1 - 1 | ||||||
| One prior stressor | Exposed | 23 | 157 | 12.8 | 180 | 87.2 |
| 2 - 0 | ||||||
| Two prior stressors | Unexposed | 30 | 60 | 33.3 | 90 | 66.7 |
| 2 - 1 | ||||||
| Two prior stressors | Exposed | 20 | 70 | 22.2 | 90 | 77.8 |
| Chi-Square Test Results - History × Pathogen Interaction | |
| Test for independence between history and pathogen exposure | |
| Statistic | Value |
|---|---|
| Chi-square statistic | 18.392800 |
| Degrees of freedom | 2.000000 |
| P-value | 0.000101 |
| Logistic Regression Results - History × Pathogen Interaction | |||||
| Model: Outcome ~ History * Pathogen | |||||
| Term | Estimate | Std_Error | z_value | p_value | p_signif |
|---|---|---|---|---|---|
| (Intercept) | 26.566 | 39,569.344 | 0.001 | 0.999 | ns |
| History1 | 0.000 | 49,279.608 | 0.000 | 1.000 | ns |
| History2 | 0.000 | 60,658.574 | 0.000 | 1.000 | ns |
| Pathogen1 | 0.000 | 55,788.569 | 0.000 | 1.000 | ns |
| History1:Pathogen1 | 0.000 | 69,158.557 | 0.000 | 1.000 | ns |
| History2:Pathogen1 | 0.000 | 83,891.969 | 0.000 | 1.000 | ns |
| Effect Size Analysis - History × Pathogen Interaction | |||||
| Odds ratios comparing exposed vs unexposed within each history level | |||||
| Comparison | Odds_Ratio | Mortality_Exposed | Mortality_Unexposed | Mortality_Difference | Interpretation |
|---|---|---|---|---|---|
| History 0: Exposed vs Unexposed | 0.9 | 8.9 | 10.0 | −1.1 | Lower mortality when exposed |
| History 1: Exposed vs Unexposed | 0.7 | 12.8 | 18.3 | −5.5 | Lower mortality when exposed |
| History 2: Exposed vs Unexposed | 0.6 | 22.2 | 33.3 | −11.1 | Lower mortality when exposed |
| Trend Analysis - Pathogen Effect Across History Levels | |
| Testing whether the effect of pathogen exposure changes across history levels | |
| Statistic | Value |
|---|---|
| Correlation coefficient (r) | -0.9976 |
| P-value | 0.044036 |
| Interpretation | Significant trend in pathogen effect across history levels |
| Pairwise Chi-Square Tests - Within Unexposed Group | |||
| Comparing history levels within unexposed fish | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| History 0 vs History 2 (Unexposed) | 13.0933 | 0.000296 | Yes |
| History 1 vs History 2 (Unexposed) | 6.7314 | 0.009473 | Yes |
| History 0 vs History 1 (Unexposed) | 2.5693 | 0.108955 | No |
| Pairwise Chi-Square Tests - Within Exposed Group | |||
| Comparing history levels within exposed fish | |||
| Comparison | Chi_square | P_value | Significant |
|---|---|---|---|
| History 0 vs History 2 (Exposed) | 5.1175 | 0.023686 | Yes |
| History 1 vs History 2 (Exposed) | 3.3228 | 0.068326 | No |
| History 0 vs History 1 (Exposed) | 0.5512 | 0.457833 | No |
| Pathogen Effect Within Each History Level | ||||||
| Comparing percent mortality between unexposed vs exposed within the same history of stressors | ||||||
| Comparison | Unexposed_Mortality | Exposed_Mortality | Mortality_Difference | Z_statistic | P_value | Significant |
|---|---|---|---|---|---|---|
| History 2 (Unexposed) vs History 2 (Exposed) | 33.3 | 22.2 | −11.1 | −1.6759 | 0.0938 | No |
| History 1 (Unexposed) vs History 1 (Exposed) | 18.3 | 12.8 | −5.5 | −1.4440 | 0.1487 | No |
| History 0 (Unexposed) vs History 0 (Exposed) | 10.0 | 8.9 | −1.1 | −0.2523 | 0.8008 | No |